﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Numerics;
using System.Security.Cryptography;
using System.Text;

namespace MainPaper
{
    internal static class AlgoAndMath
    {
        public static ulong DotProduct(ulong[] a, ulong[] b, ulong modulo)
        {
            if (a.Length != b.Length) throw new ArgumentException();
            ulong res = 0;
            for (int i = 0; i < a.Length; i++)
            {
                res = (res + a[i] * b[i]) % modulo;
            }
            return res;
        }

        public static BigInteger DotProduct(BigInteger[] a, BigInteger[] b, BigInteger modulo)
        {
            if (a.Length != b.Length) throw new ArgumentException();
            BigInteger res = 0;
            for (int i = 0; i < a.Length; i++)
            {
                res = (res + a[i] * b[i]) % modulo;
            }
            return res;
        }

        public static ulong SquareAndMultiply(ulong baseValue, ulong power, ulong modulo)
        {
            ulong res = 1;
            for (int i = sizeof (ulong)*8 - 1; i >= 0; i--)
            {
                res = FastMultilication(res, res, modulo);
                if ((power & (1UL << i)) != 0)
                {
                    res = FastMultilication(res, baseValue, modulo);
                }
            }
            return res;
        }

        public static void SplitInteger(ulong value, out ulong first, out ulong second)
        {
            first = second = 0UL;
            const int halfSize = sizeof (ulong)*4;
            for (int i = 0; i < halfSize; i++) //  sizeof(ulong) * 8 / 2
            {
                if ((value & (1UL << i)) != 0) second |= 1UL << i;
                if ((value & (1UL << (i + halfSize))) != 0) first |= 1UL << i;
            }
        }

        internal static void SplitBigInteger(BigInteger value, out BigInteger first, out BigInteger second)
        {
            var valueArr = value.ToByteArray();
            var firstArr = new byte[valueArr.Length / 2 + 1];
            var secondArr = new byte[valueArr.Length - valueArr.Length/2 + 1];
            Array.Copy(valueArr, 0, firstArr, 0, valueArr.Length / 2);
            Array.Copy(valueArr, valueArr.Length / 2, secondArr, 0, valueArr.Length - valueArr.Length / 2);
            first = new BigInteger(firstArr);
            second = new BigInteger(secondArr);
        }

        public static ulong FastMultilication(ulong a, ulong b, ulong modulo)
        {
            // replace with smth more clever!
            return (a*b)%modulo;
        }

        public static BigInteger RandomBigInt(Random generator, BigInteger upperBound)
        {
            BigInteger res;
            var buffer = upperBound.ToByteArray();
            var n = buffer.Length; // max # of bytes
            if (buffer[n - 1] == 0) n--; // remove leading zero byte
            var randomValueLength = generator.Next(n) + 1; // # of bytes in the generated value 1..n
            var highestByteValue = buffer[n - 1];

            for (int i = 0; i < buffer.Length; i++)
            {
                if (i < randomValueLength)
                {
                    buffer[i] = (byte) generator.Next(256);
                }
                else buffer[i] = 0;
            }
            if (randomValueLength == n) buffer[n - 1] = (byte) generator.Next(highestByteValue);
            res = new BigInteger(buffer);
            return res;
        }

        //http://rosettacode.org/wiki/Miller-Rabin_primality_test#C.23
        public static bool IsProbablePrime(BigInteger source, int certainty)
        {
            if (source == 2 || source == 3)
                return true;
            if (source < 2 || source%2 == 0)
                return false;

            BigInteger d = source - 1;
            int s = 0;

            while (d%2 == 0)
            {
                d /= 2;
                s += 1;
            }

            // There is no built-in method for generating random BigInteger values.
            // Instead, random BigIntegers are constructed from randomly generated
            // byte arrays of the same length as the source.
            RandomNumberGenerator rng = RandomNumberGenerator.Create();
            byte[] bytes = new byte[source.ToByteArray().LongLength];
            BigInteger a;

            for (int i = 0; i < certainty; i++)
            {
                do
                {
                    // This may raise an exception in Mono 2.10.8 and earlier.
                    // http://bugzilla.xamarin.com/show_bug.cgi?id=2761
                    rng.GetBytes(bytes);
                    a = new BigInteger(bytes);
                } while (a < 2 || a >= source - 2);

                BigInteger x = BigInteger.ModPow(a, d, source);
                if (x == 1 || x == source - 1)
                    continue;

                for (int r = 1; r < s; r++)
                {
                    x = BigInteger.ModPow(x, 2, source);
                    if (x == 1)
                        return false;
                    if (x == source - 1)
                        break;
                }

                if (x != source - 1)
                    return false;
            }

            return true;
        }

        public static BigInteger NextPrime(BigInteger p)
        {
            if (p%2 == 0) p += 1;
            else p += 2;
            for (; !IsProbablePrime(p, 40); p += 2)
            {
            }
            return p;
        }

        public static BigInteger NextSafePrime(BigInteger p)
        {
            do
            {
                p = NextPrime(p);
            } while (!IsProbablePrime((p - 1) / 2, 40));
            return p;
        }

        // Modular multiplicative inverse using extended Euclidean algorithm.
        public static BigInteger ModInverse(BigInteger a, BigInteger b)
        {
            BigInteger dividend = a % b;
            BigInteger divisor = b;

            BigInteger last_x = BigInteger.One;
            BigInteger curr_x = BigInteger.Zero;

            while (divisor.Sign > 0)
            {
                BigInteger quotient = dividend / divisor;
                BigInteger remainder = dividend % divisor;
                if (remainder.Sign <= 0)
                {
                    break;
                }

                /* This is quite clever, in the algorithm in form
                 * ax + by = gcd(a, b) we only keep track of the
                 * value curr_x and the last_x from last iteration,
                 * the y value is ignored anyway. After remainder
                 * runs to zero, we get our inverse from curr_x */
                BigInteger next_x = last_x - curr_x * quotient;
                last_x = curr_x;
                curr_x = next_x;

                dividend = divisor;
                divisor = remainder;
            }

            if (divisor != BigInteger.One)
            {
                throw new Exception("Numbers a and b are not relatively primes");
            }
            return (curr_x.Sign < 0 ? curr_x + b : curr_x);
        }
        
    }

    class LagrangeInterpolation
    {
        private readonly Point<BigInteger>[] _points;


        public LagrangeInterpolation(Point<BigInteger>[] points)
        {
            _points = points;
        }

        public BigInteger ComputeAtPoint(BigInteger point, BigInteger modulo)
        {
            var res = BigInteger.Zero;
            for (int i = 0; i < _points.Length; i++)
            {
                var top = BigInteger.One;
                var bottom = BigInteger.One;
                for (int j = 0; j < _points.Length; j++)
                {
                    if (i == j) continue;
                    var t = (point - _points[j].X);
                    if (t < 0) t += modulo;
                    top *= t;
                    top %= modulo;
                    t = (_points[i].X - _points[j].X);
                    if (t < 0) t += modulo;
                    bottom *= t;
                    bottom %= modulo;
                }
                //top %= modulo;
                //bottom %= modulo;
                //var gcd = BigInteger.GreatestCommonDivisor(top, bottom);
                //if (gcd != 1)
                //{
                //    top /= gcd;
                //    bottom /= gcd;
                //}
                bottom = AlgoAndMath.ModInverse(bottom, modulo);
                res += _points[i].Y * top * bottom;
            }
            res = res % modulo;
            if (res < 0) res += modulo;
            return res;
        }
    }

    struct Point<T>
    {
        private readonly T _x;
        private readonly T _y;

        public T X
        {
            get { return _x; }
        }

        public T Y
        {
            get { return _y; }
        }

        public Point(T x, T y)
        {
            _x = x;
            _y = y;
        }
    }
}
